library(terra)terra 1.7.71
The rgrass package allow us to interact with GRASS GIS tools (and data) serving as an interface between GRASS GIS and R. The rgrass package is developed and maintained by Roger Bivand and can be found at: https://github.com/rsbivand/rgrass/. In this fast track tutorial, we will learn how to use GRASS GIS from R.
Setup note:
To run this tutorial locally you should have GRASS GIS 8.4+, R and, optionally, RStudio installed. You will also need to install rgrass, terra and mapview R packages and download the North Carolina sample dataset.
The main functions within rgrass are the following:
initGRASS(): starts a GRASS GIS session from R.execGRASS(): executes GRASS GIS commands from R.gmeta(): prints GRASS GIS session metadata like database, project, mapset, computational region settings and CRS.read_VECT() and read_RAST(): read vector and raster maps from a GRASS project into R.write_VECT() and write_RAST(): write vector and raster objects from R into a GRASS GIS project.For further details on rgrass functionality, usage examples and data format coercion, see: https://rsbivand.github.io/rgrass/.
If you are a regular R user that needs to use GRASS GIS functionality because, well, you know it rocks, rgrass has your back. For example, maybe you struggle with large raster datasets in R or you need some specific tool, like watershed delineation for a large high resolution DEM. We will show here the way to use GRASS GIS tools within your R workflows.
On the other hand, if you already use GRASS as your geospatial data processing engine, you most likely have your spatial data within GRASS projects. You might need however to do some statistical analysis, some modelling and prediction or create publication ready visualizations in R. In such cases, you can start a GRASS GIS session in your project from R or RStudio.
Let’s see the general basic steps and then dive into the details:
rgrass library with library(rgrass)initGRASS()execGRASS()read_VECT(), read_RAST(), write_VECT() and write_RAST() to read data from and write data into GRASS database.GRASS raster and vector maps are translated into terra’s packege SpatRaster and SpatVector objects, respectively. These objects can then, within R, be easily coerced to other types of spatial objects such as simple features (sf), stars, etc.
See terra vignettes with further explanations and examples: https://rspatial.github.io/terra/.
In the case you need to include some of the cool GRASS GIS tools within your R workflow, the initGRASS() function allows you to create temporary GRASS projects to use GRASS modules on R objects. This is equivalent to what QGIS does when you use GRASS tools via the Processing Toolbox.
So, the workflow goes like this: we will use initGRASS() to create a temporary GRASS project based on the extent, resolution and CRS of a raster or vector R object, likely the one we want to process or one that has the extent of our study area. Hence, we need to pass a reference spatial grid via the SG parameter. Then, we’ll write our R objects into the temporary GRASS project, run the desired processes, export the outputs back to the R environment and done.
Let’s start with getting some spatial data, eg. a raster file, into R.
library(terra)terra 1.7.71
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
plot(r)
Now, we’ll load the rgrass library and start GRASS GIS with SG = r, so a GRASS project is internally created with r’s object CRS (BTW, you can check that with crs(r)), extent and resolution. These latter define the GRASS GIS computational region that will affect all raster processing, i.e., all new raster maps generated within GRASS GIS will have the same extent and resolution of the map provided. If you wish to change the computational region later on, you can use the g.region GRASS tool with execGRASS("g.region --h").
library(rgrass)GRASS GIS interface loaded with GRASS version: (GRASS not running)
Note that we also pass tempdir() as the folder where our GRASS temporary project will be created. tempdir() will default to whichever is the temporary folder in your operative system.
# path to GRASS binaries
grassbin <- system("grass --config path", intern = TRUE)
# start GRASS GIS from R
initGRASS(gisBase = grassbin,
home = tempdir(),
SG = r,
override = TRUE)gisdbase /tmp/Rtmp7Wq6Wd
location file5cda8663a8f01
mapset file5cda8325d3f4f
rows 90
columns 95
north 50.19167
south 49.44167
west 5.741667
east 6.533333
nsres 0.008333333
ewres 0.008333326
projection:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Now, we can write our SpatRaster into the GRASS GIS temporary project.
write_RAST(r, "terra_elev")Importing raster map <terra_elev>...
0% 3% 6% 10% 13% 16% 20% 23% 26% 30% 33% 36% 40% 43% 46% 50% 53% 56% 60% 63% 66% 70% 73% 76% 80% 83% 86% 90% 93% 96% 100%
SpatRaster read into GRASS using r.in.gdal from file
Alternatively, we can use GRASS GIS importing tools to import common raster and vector formats. These will provide re-projection on the fly, if needed.
execGRASS("r.import", input=f, output="test")Importing raster map <test>...
0% 3% 6% 10% 13% 16% 20% 23% 26% 30% 33% 36% 40% 43% 46% 50% 53% 56% 60% 63% 66% 70% 73% 76% 80% 83% 86% 90% 93% 96% 100%
Let’s check both are indeed within the project and run the GRASS module r.slope.aspect on it.
execGRASS("g.list", type = "raster")terra_elev
test
execGRASS("r.slope.aspect",
elevation = "terra_elev",
slope = "slope",
aspect = "aspect") 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 34% 37% 40% 43% 46% 49% 52% 56% 59% 62% 65% 68% 71% 74% 78% 81% 84% 87% 90% 93% 96% 100%
Aspect raster map <aspect> complete
Slope raster map <slope> complete
execGRASS("g.list", type = "raster")aspect
slope
terra_elev
test
Let’s get slope and aspect maps into R
grass_maps <- read_RAST(c("aspect", "slope"))Checking GDAL data type and nodata value...
2% 5% 8% 11% 14% 17% 20% 23% 26% 30% 33% 36% 40% 43% 46% 50% 53% 56% 60% 63% 66% 70% 73% 76% 80% 83% 86% 90% 93% 96% 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
2% 5% 8% 11% 14% 17% 20% 23% 26% 30% 33% 36% 40% 43% 46% 50% 53% 56% 60% 63% 66% 70% 73% 76% 80% 83% 86% 90% 93% 96% 100%
r.out.gdal complete. File </tmp/Rtmp7Wq6Wd/file5cda866f11586.grd> created.
Checking GDAL data type and nodata value...
2% 5% 8% 11% 14% 17% 20% 23% 26% 30% 33% 36% 40% 43% 46% 50% 53% 56% 60% 63% 66% 70% 73% 76% 80% 83% 86% 90% 93% 96% 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
2% 5% 8% 11% 14% 17% 20% 23% 26% 30% 33% 36% 40% 43% 46% 50% 53% 56% 60% 63% 66% 70% 73% 76% 80% 83% 86% 90% 93% 96% 100%
r.out.gdal complete. File </tmp/Rtmp7Wq6Wd/file5cda858f9b728.grd> created.
grass_mapsclass : SpatRaster
dimensions : 90, 95, 2 (nrow, ncol, nlyr)
resolution : 0.008333326, 0.008333333 (x, y)
extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
sources : file5cda866f11586.grd
file5cda858f9b728.grd
names : aspect, slope
min values : 0.08974174, 0.01416342
max values : 360.00000000, 7.22943783
Now that the output maps are back into our R environment, we can plot them, do further analysis or write them into other raster formats, in which case we use terra::writeRaster() function.
plot(grass_maps)
writeRaster(grass_maps, "grass_maps.tif", overwrite=TRUE)Alternatively, we can use GRASS GIS exporting tools like r.out.gdal and v.out.ogr, to directly save our outputs into common raster or vector files respectively.
execGRASS("r.out.gdal", input="slope", output="slope.tif", format="GTiff", flags="overwrite")Checking GDAL data type and nodata value...
2% 5% 8% 11% 14% 17% 20% 23% 26% 30% 33% 36% 40% 43% 46% 50% 53% 56% 60% 63% 66% 70% 73% 76% 80% 83% 86% 90% 93% 96% 100%
Using GDAL data type <Float32>
Input raster map contains cells with NULL-value (no-data). The value nan
will be used to represent no-data values in the input map. You can specify
a nodata value with the nodata option.
Exporting raster data to GTiff format...
ERROR 6: slope.tif, band 1: SetColorTable() only supported for Byte or UInt16 bands in TIFF format.
2% 5% 8% 11% 14% 17% 20% 23% 26% 30% 33% 36% 40% 43% 46% 50% 53% 56% 60% 63% 66% 70% 73% 76% 80% 83% 86% 90% 93% 96% 100%
r.out.gdal complete. File <slope.tif> created.
Let’s see an example for the case when we do all our geospatial data processing within GRASS and hence have all the spatial data organized within GRASS projects but we need to run some statistical analysis, modelling, prediction or visualization in R.
We start R or Rstudio and load the rgrass library. It will tell us that GRASS is not running, but we know that already… and that’s about to change in a moment.
library(rgrass)So, we start GRASS from within R or RStudio using the initGRASS() function. Since we want to start GRASS in a specific project and mapset, we need to specify them. Optionally, we specify which GRASS binary to use. This might be useful in the case we have many GRASS versions in our system. If not provided, initGRASS() will attempt to find it in default locations depending on your operative system.
We can now list and read our GRASS raster and vector maps into R and do our statistical analysis, modelling and/or visualizations using other R packages. Obviously, we can also use rgrass as a mere interface to analyze data within GRASS, i.e., no reading from or writing to GRASS GIS needed, we can just use GRASS from R. Here, we’ll demonstrate the use of all the main rgrass functions mentioned above.
Let’s then list our GRASS raster and vector maps:
# list GRASS raster maps
execGRASS("g.list", parameters = list(type="raster"))basins
developed
elevation
elevation_2
elevation_shade
geology
lakes
landuse
result_from_R
soils
# list GRASS vector maps
execGRASS("g.list", parameters = list(type="vector"))boundary_region
boundary_state
census
elev_points
firestations
geology
geonames
hospitals
points_of_interest
railroads
roadsmajor
schools
streams
streets
zipcodes
The resulting map lists could be saved in an R object that we can subset later in case we want to import several but not all raster maps, for example. Let’s see how to do that.
# save map list in an object
rast_list <- execGRASS("g.list", parameters = list(type="raster"))basins
developed
elevation
elevation_2
elevation_shade
geology
lakes
landuse
result_from_R
soils
rast_list[1] 0
attr(,"resOut")
[1] "basins" "developed" "elevation" "elevation_2"
[5] "elevation_shade" "geology" "lakes" "landuse"
[9] "result_from_R" "soils"
attr(,"resErr")
character(0)
# retrieve only the maps list and overwrite rast_list
rast_list <- attributes(rast_list)$resOut
# import elevation and landuse
to_import <- c("elevation", "landuse") # optionally, by position: rast_list[c(3,7)]
maplist <- list()
for (i in to_import) {
maplist[i] <- read_RAST(i)
}Checking GDAL data type and nodata value...
2% 5% 8% 11% 14% 17% 20% 23% 26% 29% 32% 35% 38% 41% 44% 47% 50% 53% 56% 59% 62% 65% 68% 71% 74% 77% 80% 83% 86% 89% 92% 95% 98% 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
2% 5% 8% 11% 14% 17% 20% 23% 26% 29% 32% 35% 38% 41% 44% 47% 50% 53% 56% 59% 62% 65% 68% 71% 74% 77% 80% 83% 86% 89% 92% 95% 98% 100%
r.out.gdal complete. File </tmp/Rtmp7Wq6Wd/file5cda8613bd54d.grd> created.
Warning in `[<-`(`*tmp*`, i, value = read_RAST(i)): implicit list embedding of
S4 objects is deprecated
Checking GDAL data type and nodata value...
2% 5% 8% 11% 14% 17% 20% 23% 26% 29% 32% 35% 38% 41% 44% 47% 50% 53% 56% 59% 62% 65% 68% 71% 74% 77% 80% 83% 86% 89% 92% 95% 98% 100%
Using GDAL data type <UInt32>
Exporting raster data to RRASTER format...
2% 5% 8% 11% 14% 17% 20% 23% 26% 29% 32% 35% 38% 41% 44% 47% 50% 53% 56% 59% 62% 65% 68% 71% 74% 77% 80% 83% 86% 89% 92% 95% 98% 100%
r.out.gdal complete. File </tmp/Rtmp7Wq6Wd/file5cda8a157673.grd> created.
Warning in `[<-`(`*tmp*`, i, value = read_RAST(i)): implicit list embedding of
S4 objects is deprecated
maplist$elevation
class : SpatRaster
dimensions : 1350, 1500, 1 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
extent : 630000, 645000, 215000, 228500 (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(HARN) / North Carolina (EPSG:3358)
source : file5cda8613bd54d.grd
name : file5cda8613bd54d
min value : 55.57879
max value : 156.32986
$landuse
class : SpatRaster
dimensions : 1350, 1500, 1 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
extent : 630000, 645000, 215000, 228500 (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(HARN) / North Carolina (EPSG:3358)
source : file5cda8a157673.grd
name : file5cda8a157673
min value : 0
max value : 7
Remember that raster objects will always be exported from GRASS GIS following the computational region settings. So, bear that in mind when reading into R which will hold them in memory. Vectors however will be exported in their full extent.
Let’s load the terra library to quickly display our recently imported raster maps:
library(terra)
plot(maplist$elevation)
Optionally, we could stack our two SpatRaster objects together and plot them together:
rstack <- rast(maplist)
plot(rstack)
Let’s create a boxplot of elevation per land class.
boxplot(rstack$elevation, rstack$landuse, maxcell=50000)
Let’s import a vector map, too, and explore its attributes.
census <- read_VECT("census")Exporting 2547 areas (may take some time)...
5% 11% 17% 23% 29% 35% 41% 47% 53% 59% 65% 71% 77% 83% 89% 95% 100%
WARNING: 10 features without category were skipped. Features without
category are written only when -c flag is given.
v.out.ogr complete. 2537 features (Polygon type) written to <census> (GPKG
format).
head(census) cat OBJECTID AREA PERIMETER BLOCK_ BLOCK_ID BLOCKNUM TOTAL_POP
1 86 82802 9981.633 444.547 82803 82802 371830535104051 13
2 88 82805 304.076 91.371 82806 82805 371830535104046 0
3 91 82811 2400.879 192.304 82812 82811 371830535104034 20
4 140 82955 1239.619 147.011 82956 82955 371830535104039 0
5 173 83073 10358.715 429.738 83074 83073 371830535014007 80
6 228 83352 15459.671 458.961 83353 83352 371830535015007 18
POP_1RACE WHITE_ONLY BLACK_ONLY AMIND_ONLY ASIAN_ONLY HWPAC_ONLY OTHER_ONLY
1 13 13 0 0 0 0 0
2 0 0 0 0 0 0 0
3 20 3 0 0 17 0 0
4 0 0 0 0 0 0 0
5 76 38 9 0 0 0 29
6 18 17 1 0 0 0 0
POP_2RACES HISPANIC MALE FEMALE P_UNDER_5 P5_TO_9 P10_TO_14 P15_TO_19
1 0 0 6 7 1 2 1 1
2 0 0 0 0 0 0 0 0
3 0 0 10 10 1 1 1 1
4 0 0 0 0 0 0 0 0
5 4 54 37 43 20 4 5 7
6 0 3 10 8 1 0 1 0
P20_TO_24 P25_TO_34 P35_TO_44 P45_TO_54 P55_TO_59 P60_TO_64 P65_TO_69
1 1 3 1 1 1 0 1
2 0 0 0 0 0 0 0
3 2 8 2 3 0 1 0
4 0 0 0 0 0 0 0
5 14 17 7 5 0 0 0
6 0 5 0 5 1 2 2
P70_TO_74 P75_TO_84 P75_OLDER P85_OLDER MEDIAN_AGE P18_OLDER P21_OLDER
1 0 0 0 0 25.5 9 8
2 0 0 0 0 0.0 0 0
3 0 0 0 0 30.0 16 16
4 0 0 0 0 0.0 0 0
5 0 1 1 0 21.8 48 44
6 0 1 1 0 52.5 16 16
P65_OLDER M75_OLDER F75_OLDER HOUSEHOLDS HH_SIZE SINGLE_HH FAMILY_HH
1 1 0 0 3 4.33 0 2
2 0 0 0 0 0.00 0 0
3 0 0 0 7 2.86 0 6
4 0 0 0 0 0.00 0 0
5 1 0 1 17 4.71 1 16
6 3 0 1 9 2.00 3 6
MAR_COUPLE W_OWN_CHLD M_SNGL_PAR F_SNGL_PAR NONFAM_HH HH_W_CHILD HH_W_ELDER
1 2 2 0 0 1 3 1
2 0 0 0 0 0 0 0
3 6 3 0 0 1 3 0
4 0 0 0 0 0 0 0
5 5 3 3 3 0 10 1
6 6 2 0 0 0 2 2
SNGL_ELDER POP_IN_HH FAM_HHLDER SPOUSE CHILD GROUP_QTRS GQ_INST GQ_NONINST
1 0 13 2 2 4 0 0 0
2 0 0 0 0 0 0 0 0
3 0 20 6 6 6 0 0 0
4 0 0 0 0 0 0 0 0
5 0 80 16 5 22 0 0 0
6 1 18 6 6 3 0 0 0
FAM_SIZE HOUSING_U OCCUPIED VACANT OWNER_U RENTER_U VACANT_SEA HH_SIZE_OW
1 4.00 3 3 0 0 3 0 0.00
2 0.00 0 0 0 0 0 0 0.00
3 3.00 9 7 2 7 0 0 2.86
4 0.00 0 0 0 0 0 0 0.00
5 3.81 17 17 0 9 8 0 5.00
6 2.50 9 9 0 8 1 0 1.88
HH_SIZE_RE
1 4
2 0
3 0
4 0
5 4
6 3
summary(census$TOTAL_POP) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 0.00 26.00 63.03 62.00 6173.00
plot(census, "P25_TO_34", type="interval", breaks=5, plg=list(x="topright"))
Let’s do some interactive visualization with mapview.
library(mapview)
mapview(rstack$elevation) + census